Bulletin of Mathematical Biology manuscript No. 

(will be inserted by the editor) 



Maximal Sensitive Dependence and the Optimal Path to 
Epidemic Extinction 

Eric Forgoston • Simone Bianco • Leah B. 
Shaw ■ Ira B. Schwartz 



Received: date / Accepted: date 



Abstract Extinction of an epidemic or a species is a rare event that occurs due to a 
large, rare stochastic fluctuation. Although the extinction process is dynamically un- 
stable, it follows an optimal path that maximizes the probability of extinction. We 
show that the optimal path is also directly related to the finite-time Lyapunov expo- 
nents of the underlying dynamical system in that the optimal path displays maximum 
sensitivity to initial conditions. We consider several stochastic epidemic models, and 
examine the extinction process in a dynamical systems framework. Using the dynamics 
of the finite-time Lyapunov exponents as a constructive tool, we demonstrate that the 
dynamical systems viewpoint of extinction evolves naturally toward the optimal path. 

Keywords Stochastic dynamical systems and Lyapunov exponents ■ Optimal path 
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1 Introduction 

Control and eradication of infectious diseases are among the most important goals for 
improving public health. Although the global eradication of a disease (e.g. smallpox) 
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has been achieved ( Breman and Arital . ll98oh . it is difficult to accomplish and remains 
a public health target for polio , malaria, and many other diseases, including childhood 
diseases ( Avlward et all l2000h . The global eradication of a disease is not the only type 
of disease extinction process. For example, a disease may "fade out" or become locally 
extinct in a re gion. Since the extin ction is local, the disease may be reintroduced from 
other regions (jCrasslv et alll2005h . Additionally, extinction may also occur to individ- 
ual strains of a multistrain disease ( Minavev and Fergusonl 120091 ) , such as influenza or 
dengue fever. It is worth noting that the extinction process is also of interest in the fields 
of evolution and ecology. As an example, bio -diversity arises from the interplay between 
the i ntroduction and extinction of species ( Azaele et all 120061 : iBanavar and Maritanl . 
120091 ). 

In order to predict the dynamics of disease outbreak and spread as well as to 
implement control strategies which promote disease extinction, one must investigate 
how m odel parameters affect the probability of extinction. For example. iDvkman et al.l 
( 2008j) showed that disease control and extinction depend on both epidemiological and 
sociological parameters determining epidemic growth rate. Additionally, since extinc- 
tion occurs in finite populations, another factor in determining extinction risk is the lo- 
cal co mmunity size (|Bartlettl . E957l . ll960l : lkeeling and Grenfeljfl997l : IConlan and Grenfeil . 
l2007h . Further complications arise from the fact that in general, stochastic effects cause 
random transitions within the discrete, finite populations. These stochastic effects may 
be intrinsic to the system or may arise from the external environment. Small popula- 
tion size, low contact frequency for fre quency-dependent transmiss ion, competition for 
resources, and evolutionary pressure Jde Castro and Bolkeri l2005h . as well as hetero- 
geneity in populations and transmission (|LloTd~eTldTT2*0*0*7l ) may all be determining 
factors for extinction to occur. Other factors which can affect the risk of e xtinction 
include the nature a nd strength o f the n oise ( Melbourne and Hastings! 120081') , disease 
outb reak amplitude ( Alonso et all 120061 ') , and seasonal phase occurrence (|Stone et all 
120071 '). 

In large populations, the intensity of the intrinsic noise is often quite small. How- 
ever, it is possible that a rare, large fluctuation which occurs with finite probability can 
cause the system to reach the extinct state. Since the extinct state is absorbing due 
to effective stocha stic forces, eventual extinction is guaranteed when t here i s no source 
of reintroduction (|Bartlettl. Il949l ; lAllen and Burginl . |2000| ; iGardinerl . |2004 ). However, 
because fade outs are usually rare events in large populations, typical time scales for 
extinction may be extre mely long. 

Birth-death systems (|Gardinerj , |2004| ; Ivan Kampenl, 120071 ) , which possess intrinsic 
noise, form an important class of stochas tic processes. These systems have been used in 
the field of mathem atical epidemiology ( Bartlettl . Il96ll : lAndersson and Brittonl . l2000l : 
IChoisv et afll2007l ). In these systems, the intrinsic noise arises from the discreteness of 
the individuals (or species) and the randomness of their interactions. To predict prob- 
abilities of events occurring at certain times, a description of the stochastic system is 
provided using the master equation formalism. Stochastic master equations ar e com- 
monl y used in statistical physics when dealing with chemical reaction processes (|Kubol . 
119631 1. 

For systems with a large number of individuals, a van Kampen system-size ex- 
pansi on may be used to approximate th e master equation by a Fokker-Planck equa- 
tion ( Gardinerl . 120041 : Ivan Kampenl . |2007 | ) . However, the t echnique fails in determinin g 
the probability of large fluctuations ( Gaveau et all 1 19961 ; lElgart and Kamenevl 120041 ') . 
Instead, in this article, we will employ an eikonal approximation to recast the problem 
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in terms of an effective classical Hamiltonian system ( Kamenev and Meersonl . I2008T 1. 
With high probability, the intrinsic noise will induce extinction of the disease or species 
along a heteroclinic trajectory in the phase space of the classical Hamiltonian flow. This 
trajectory is known as the optimal path (to extinction). 

It is highly desirable to locate the optimal path since the extinction rate depends 
on the probability to traverse this path. Additionally, the effect of a co ntrol strategy on 
the e xtinction rate can be determined by its effect on the optimal path ( Dvkman et al.l . 
2008J). Through the use of the eikonal approximation, we consider a mechanistic dy- 
namical systems problem rather than the original stochastic problem. Unlike other 
methods, this approach enables one both to estimate extinction lifetimes and to draw 
conclusions about the path taken to extinction. This more detailed u nderstanding of 
how extinction occurs may lead to new stochastic control strategies ( Dvkman et all 
120081 ). 

In general, the optimal path to extinction is an unstable dynamical object. There- 
fore, many researchers have investigated the scaling of extinction r ates with respect to a 
parameter near a bifurcation point, where the dynamics are slow jDoering et alj . l2005l : 
iKamenev and Meersonl . bOOsfeamenev et all l20Qg| ; iDvkman et al.!l200Sl) . The analvt- 

ical calculation of extinction rates far from a bifurcation point is much mo re difficult, 

and th us, researchers often resort to averaging over many stochastic runs fe.g. lShaw and Schwartj 

The numerical computation of the optimal path trajector y has been achieved 
in the past using a shooting method ( Kamenev and Meersonl . 120081 1. However, since the 
procedure is very sensitive to boundary conditions, it is difficult to implement when 
analyzing paths far away from bifurcation points or when analyzing high-dimensional 
models. 

In this article we develop a novel method for computing the optimal extinction 
trajectory that relies on the calculation of the dynamical system's finite-time Lyapunov 
exponents (FTLE). The classical Lyapunov exponent provides a quantitative measure 
of how infinitesim ally close particles in a dynamic al system behave asymptotically 
as time t — > ±oo ( Guckenheimer and Holmes! . Il986f) . Similarly, the FTLE provides a 
quantitative measure of how much nearby particles separate after a specific amount of 
time has elapsed. 

Th e FTLE fields were used by IPierrehumbertl ( 199ll ) and lPierrehumbert and Yand 
( 19931 ) to characterize structures, including mixing regions and transport barriers, 
in the atmosphere. Later, these structures were quantified more r i gorou s ly by intro- 
ducing the ide a of Lagrangian Coherent Structures (XCS) llHallerl. iboocl . l200ll . 120021 : 
IShadden et~al] . 120051 ; iLekien et ail . |2007l : iBranicki and Wiggins] . l2009|) . The definition 
of a LCS as a ridg e of th e FTLE field was introduced bv lHallerl ( 20021 ) and formalized 
by IShadden et all (|2005l ). The FTLE method provides a measure of how sensitively 
the system's future behavior depends on its current state, and the LCS, or ridge, is 
a structure which has locally maximal FTLE value. In this article, we will show that 
the system displays maximum sensitivity near the optimal extinction trajectory, which 
enables us to dynamically evolve toward the optimal escape trajectory using FTLE 
calculations. 

In this article, we consider three standard epidemic models that contain intrinsic or 
external noise sources and illustrate the power of our method to locate the optimal path 
to extinction. Even though our examples are taken from infectious disease models, the 
approach is applicable to any extinction process or escape process. Section [2] discusses 
stochastic modeling in the limit of large population size, while Sec. [3] describes the 
theory that underlies the Lyapunov exponent computations. In Sec. [4j our method is 
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used to find the optimal path to extinction for three examples. In Sec. [5] we demonstrate 
that the optimal path to extinction also possesses a local maximum of the FTLE field, 
and conclusions are contained in Sec. [S] 



2 Stochastic modeling in the large population limit 

To introduce the idea of constructing an optimal path in stochastic dynamical systems, 
we consider the problem of extinction taken from mathematical epidemiology. The 
stochastic formulation of the problem accounts for the random state transitions which 
govern the dynamics of the system. To quantitatively account for the randomness in 
the system, we will formulate a master equation which describes the time e volution of 
the probability of f i nding the system in a particular state at a certain time ( Gardiner! . 
120041 : Ivan KampenL 120071 ). 

The state variables X £ R™ of the system describe the components of a popu- 
lation, while the random state transitions which govern the dynamics are described 
by the transition rates W(X,r), where r £ ffi™ is an increment in the change of X. 
Consideration of the net change in increments of the state of the system results in the 
following master equation for the probability density p(X ,t) of finding the system in 
state X at time t: 

= £ [W(X - r; r)p(X - r, t) - W(X;r)p(X, t)) . (1) 

r 

If the total population size of the system is N, the components of the state variable 
can be rescaled to be fractions of the population by letting x — X/N. Then the 
general solution ( Kubo et al.l . Il973l ) for the transition probability from x$ to x in the 



time interval from to to t can be written as the following path integral: 

P(x,t\x ,t ) = / dV(x,p)e X p\-N / ds [H(x(s),p(s),s) - p(s)x(s)] \ , (2) 



to 



where dT>(x,p) denotes integration over all paths. 

The integral in the exponent of Eq. @ is the action, and the Hamiltonian H (x, p; t) 
is given in general as 

H(x,p-t) = ^2w(x;r)(exp(p-r) - 1), (3) 

r 

where w(x;r) is defined as the transition rate W per individual. The Hamiltonian H 
depends both on the state of the system x as well as the momentum p, which provides 
an effective force due to stochastic fluctuations on the system. It should be noted 
that instead of using the Hamiltonian representation, one could use the Lagrangian 
representation, which results in the following alternative solution: 

P(sB,t|a;o>to) = I dV(x)oxp ^ N J ds L(x(s),x(s), s)^j , (4) 



with L(x, x; t) = —H(x,p; t) + x ■ p. 
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The action reveals much information about the probability evolution of the system, 
from scaling near bifurcation points in non- Gaussian proce s ses to rates of extinc tion 
as a function of epidemiological parameters I Dvkmarl Il990l ; iDvkman et al, 1, 120081 ). In 
order to maximize the probability of extinction, one must minimize the action. The 
minimizing formulation entails finding the solution to the Hamilton- Jacobi equation, 
which means that one must solve the 2n-dimensional system of Hamilton's equations 
for x and p. The appropriate boundary conditions of the system are such that a solution 
starts at a non-zero state, such as an endemic state, and asymptotically approaches 
one or more zero components of the state vector. Therefore, a trajectory that is a 
solution to the two-point boundary value problem determines a path, which in turn 
yields the probability of going from the initial state to the final state. The optimal 
path to extinction is the path that minimizes the action in either the Hamiltonian or 
Lagrangian representation. 

To illustrate an application of the theory for a finite population, we consider in de- 
tail an explicit example, namely the well-known problem of extinction in a Susceptible- 
Infectious- Susceptible (SIS) epidemic model. In this example, the population consists of 
S susceptible individuals and I infectious individuals. The population is driven via the 
birth rate p, which is also equal to the death rate. The transition rates for X — (S, I) T 
are given as follows: 



W(X;(l,0))=Np, W(X;(-l,0))= f xX 1 , 
W(X; (0, -1)) = pX 2 , W(X; (1, -1)) = 1 X 2 , 

W(X;(-1,1)) = PX X X 2 /N, 



(5) 



where j3 is the mass action contact rate, 7 is the recovery rate, and N is now a parameter 
for the average size of the population. For large S, I oc N, fluctuations of S and I are 
small on average. If these fluctuations are disregarded, one arrives at the following 
deterministic mean-field equations for S and /: 



Xi =Nfi - fiXi + iX 2 - pXxX 2 /N, 
X 2 = -(p + -y)X 2 + pX x X 2 jN. 



(6a) 
(6b) 



Equations (|6a|) - (|6b|) are the standard equations of the SIS model in the absence of 
fluctuations. For the parameter Rq — /3/(p + 7) > 1, they have a stable, attracting 
solution Xj\ = Nxa with x\a = Rq 1 , and x 2 a = 1 — Rq 1 , which corresponds to 



endemic disease. In addition, Eqs. (|6a|) - (|6bl) have an unstable stationary state (saddle 
point) given by Xg = Nxg with x^g = 1 and x 2 g = 0, which corresponds to the 
extinct, or disease-free, state. 

For N ^> 1, the steady state distribution p(X) has a peak at the stable state 
X a with width o c N 1 ^ 2 . T his peak is formed ov er a typical relaxation time given 



Dvkman et al. I (|2008h and lSchwartz et all (|2009l ). However, in the process of extinc- 



tion, we are interested in the probability of having a small number of infectious individ- 



tained bv seeking' the solution of Eqs. (|6a|)-(16bl) in the eikonal form 


Elgart and Kamenev 


2004 


DoerinE et all 2005; Kubo et all 19731; Wentzelll.ll976!;lGane 


. 19871; Dvkman et al. 


1994 


Tretiakov et al. 20031") given bv 



p(X) = exp l-NS(x)], 

p(X + r) sa p(X) exp(-p ■ r) 



= X/N, 

p = dS(x)/dx, 



(7) 
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where S(x) is the action. 

To leading order in N~ 1 , the equation for S(x) has a form of the Hamilton-Jacobi 
equation S = —H(x,d x S;t), where S is the effective action, and the effective Hamil- 
tonian is given by Eq. ©, with w(x; r) = N~ 1 W(X;r) being the transition rates per 
individual. The action S(x) can be found from classical trajectories of the auxiliary 
system with Hamiltonian H that satisfy the following equations: 

x = d p H(x,p), p=-d x H(x,p). (8) 

Since the maximum of the probability of extinction is found by minimizing the 
action, we compute the trajectory satisfying the Hamiltonian system that has as its 
asymptotic limits in time the endemic state as t -> — oo and the extinct state as 
t -> +00. The action then has the form, f rom Eq. © dWentzell Il97rj ; iGand , Il987l ; 
iDvkman etafl . 1 19941 : iTretiakov et all 120031 ) : 



S(x s )= / p-xdt, H{x,p)=0. (9) 



In Eq. ([2}, the integral is calculated for a Hamiltonian trajectory (x(t),p(t)^ T that 
starts as t — s- —00 at x — > xj^,p — > 0, and arrives as t — > 00 at the state xg. This 
trajectory describes the most probable sequence of elementary events x — s- x + r that 
brings the system to xg. 

Several authors have consider ed how extinction ra t es scale with respect to a pa- 
rameter near b i furca t ion points dPoering et all 120051 : H amenev and Meersonl . 120081 : 
iKamenev et alll2008l : IPvkman et al. . 20081 ) when the distance to the bifurcation point 



is small and the dynamics is very slow. For an epidemic model, this means that the 
reproductive rate of infection is greater than but very close to one. However, most 
real d iseases have reproductive rates of infection greater than 1.5 ( Anderson and Mavl . 



Il99lh . which translates into faster growth rates from the extinct state. In general, in 
order to get analytic scaling results, one must compute the optimal path using either 
the Hamiltonian or Lagrangian equations of motion. However, far from bifurcation 
points, one is seldom able to perform the required analysis or computation. 

Additionally, the computation of the opt i mal p ath involves the use of a numeri- 



Additionally, tne computati on 01 tne opt i mal p atn involves tne use or a numeri- 
cal shooting method ( Kamenev and Meersonl . I2008T I. The Hamiltonian or Lagrangian 



representation of an n-dimensional dynamical system lies in 2n-dimensional space. 
Therefore, for even relatively low-dimensional dynamical systems, the use of a shoot- 
ing method to find the optimal path to extinction can be quite problematic. In the next 
section, we demonstrate how to evolve naturally to the optimal path using a dynamical 
systems approach. 



3 Finite-Time Lyapunov Exponents 

We consider a velocity field v : R 2 " x I -> R 2n given by the Hamiltonian field in Eq. (© 
which is defined over the time interval I = [t^, tf] C ffi and the following system of 
equations: 



v(t;U,Vo) = v(y{t;ti,y ),t), 
y(ti;ti,y ) = yo, 



(10a) 
(10b) 
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where y = (x,p) T € M 2n , y £ K 2 ", and t G J. 

Such a continuous time dynamical system has quantities, known as Lyapunov ex- 
ponents, which are associated with the trajectory of the system in an infinite time 
limit. The Lyapunov exponents measure the growth rates of the linearized dynam- 
ics about the trajectory. To find the finite-time Lyapunov exponents (FTLE), one 
computes the Lyapunov exponents on a restricted finite time interval. For each ini- 
tial condition, the exponents provide a measure of its sensitivity to small pertur- 
bations. Therefore, the FTLE is a measure of the local sensitivity to initial data. 
For the purpose of completeness, we briefly recapitulate the derivation of the FTLE. 
Details regardi ng the derivation along w it h the appropri a te sm o othness assumption s 



can be found in 



Hailed J2000L l200ll l2002h . IShadden etaD (|2005l l. lLekien et all (|2007f ) 



and lBranicki and Wiggins! ( 20091 ) 



The integration of Eqs. (|10a[) - (110b[l from the initial time tj to the final time t, + T 
yields the flow map (f>l' +T which is defined as follows: 

4]+ T : y h+ 4>u +T {Vo) = v(U + T; t h y ). (11) 
Then the FTLE can be defined as 



a(y, t h T) = —r In ^X max (A), (12) 

where Amax(^l) is the maximum eigenvalue of the right Cauchy-Green deformation 
tensor A, which is given as follows: 

with * denoting the adjoint. 

For a given y £ R 2n at an initial time tj, Eq. (|12|1 gives the maximum finite-time 
Lyapunov exponent for some finite integration time T (forward or backward), and 
provides a measure of the sensitivity of a trajectory to small perturbations. 

The FTLE field given by a(y,ti,T) can be shown to exhibit "ridges" of local 
maxima in phase space. The ridges of the field indicate the location of attracting 
(backward time FTLE field) and repelling (forward time FTLE field) structures. In 
two-dimensional (2D) space, the ridge is a curve which locally maximizes the FTLE 
field so that transverse to the ridge one finds the FTLE to be a local maximum. What 
is remarkable is that the FTLE ridges correspond to the optimal path trajectories, 
which is shown heuristically in Sec. [S] The basic idea is that since the optimal path 
is inherently unstable and observed only through many realizations of stochastic ex- 
periments, the FTLE shows that locally, the path is also the most sensitive to initial 
data. Figure Q] shows a schematic that demonstrates why the optimal path has a lo- 
cal maximum to sensitivity. If one chooses an initial point on either side of the path 
near the endemic state, the two trajectories will separate exponentially in time. This is 
due to the fact that both extinct and endemic states are unstable, and the connecting 
trajectory defining the path is unstable as well. Any initial points starting near the 
optimal path will leave the neighborhood in short time. 



s 




Fig. 1 (Color online) Schematic showing the path from the endemic state (blue) to the extinct 
state (red). The optimal path leaves the endemic point along an unstable manifold and connects 
to the extinct state along a stable manifold. The magenta and green dashed lines represent 
trajectories on either side of the optimal path. The initial starting distance between trajectories 
near the endemic state expands exponentially in forward time (shown by the cyan lines). 
Locally, this demonstrates that the finite-time Lyapunov measure of sensitivity with respect 
to initial data is maximal along the optimal path. This is evident in the ridges observed in the 
evolution of the exponents. 



4 Finding the Optimal Path to Extinction Using FTLE 



We now apply our theory of dynamical sensitivity to the problem of locating optimal 
paths to extinction for several examples. We consider the case of internal fluctuations, 
where the noise is not known a priori, as well as the case of external noise, where the 
noise is specified. In each case, the interaction of the noise and state of the systems 
begin with a description of the Hamiltonian (or Lagrangian) to describe the unstable 
flow. Then the corresponding equations of motion are used to compute the ridges 
corres ponding to maximum F TLE, which in turn correspond to the optimal extinction 
paths i Schwartz et alj . feoich . 



4.1 Example 1 - Extinction in a Branching- Annihilation Process 



For an example of a system with intrinsic noise fluctuations which has an analytical 
solution, we consider extinction in the stochastic branching-annihilation process 

A^2A and 2A A 0, (14) 
wher e A and fj. > are constant reaction rates l Elgart and Kamenevl . l2004l ; lAssaf et al.l 



200&f). Equation (|14|) is a single species birth-death process and can b e thought of as a 
simplified form of the Verhulst logistic model for population growth (iNaselilioOll ). 

The stochastic process given by Eq. (| 14p contains intrinsic noise which arises from 
the randomness of the reactions and the fact that the population consists of discrete 
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individuals. This intrinsic noise, which can generate a rare sequence of events that cause 
the system to evolve to the empty state, can be accounted for using a master equation 
approach. The probability P n (t) to observe, at time t, n individuals is governed by the 
following master equation: 

Pn = ~ [(n + 2)(n + l)P n+2 - n{n - l)P n ] + \[(n- l)P n -l - nP„] . (15) 



Using the formalism of lAssaf et ail (|2008h . Eq. CEHJ is recast as the following exact 
evolution equation for G(p,t): 



m=^ 1 -^W + Xp(p - 1 W (16) 



where G is a probability generating function given by 

oo 

G(p,t) = J2p np n(t), (17) 

and where p is an auxiliary variable. 

We substitute the eikonal ansatz G(p,t) = exp[— S(p, t)], where S is the action, 
into Eq. (|16|l and neglect the higher-order d 2 S/dp 2 term. This results in a Hamilton- 
Jacobi equation for S(p,t). By introducing a conjugate coordinate q = —dS/dp and 
by shifting the momentum p = p — 1, then one arrives at the following Hamiltonian: 

H{q,p)=(\{l+p)-£(2+p)q^qp. (18) 
Hamilton's equations are therefore given as: 

q=^=q[X(l + 2p)-f,(l+p)q], (19a) 

P = -^=pK2 + p)9-A(l+p)]. (19b) 

The Hamiltonian given by Eq. (|18|l has three zero-energy curves. The first is the 
mean-field zero-energy line p — 0, which contains two hyperbolic points given as hi = 
(q,p) = (A//x, 0) and ho — (q,p) = (0, 0). The second is the extinction line q = 0, which 
contains another hyperbolic point given as /12 = (q,p) = (0, —1). The third zero-energy 
curve is non-trivial and is given as follows: 

,= ^±f (20) 

The segment of the curve given by Eq. 1)201) which lies between — 1 < p < corresponds 
to a heteroclinic trajectory. As t approaches —00, the trajectory exits the hyperbolic 
point hi along its unstable manifold and enters the hyperbolic point h-2 along its 
stable manifold as t approaches 00. This heteroclinic trajectory is the optimal path to 
extinction and describes the most probable (rare) s equence of events w hich evolves the 



system from a quasi-stationary state to extinction ( Assaf et all . l20oih . 

To show that the FTLE evolves to the optimal path, we calculate the FTLE field 
using Eqs. (119a[l - ()19b[) . Figure [2ja) shows the forward FTLE plot computed using 
Eqs. (fl9li| - (fT9b)) for T = 6, with A = 2.0 and p = 0.5. In Fig. 0a), one can see that the 
optimal path to extinction is given by the ridge associated with the maximum FTLE. 
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Fig. 2 (Color online) FTLE flow fields computed using Eqs. H19al l- lll9bt with A = 2.0 and 
(j, = 0.5. The integration time is T = 6 with an integration step size of t = 0.1 and a grid 
resolution of 0.01 in both q and p. (a) Forward FTLE field with the optimal path to extinction 
given by the LCS. (b) Average of the forward and backward FTLE fields with the three zero- 
energy curves given by the LCS and overlaid with the analytical solution of these curves given 
by p = 0, q = 0, and Eq. 1201 1. Note that the averaging affects only the value of the FTLE and 
not the structure of the FTLE field. 



Not all of the attracting structures are shown in Fig. UK a) because the maximum 
forward time FTLE identifies repelling structures where nearby initial conditions di- 
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verge. The other attracting structures may be found by computing the backward FTLE 
field. By overlaying the forward and backward FTLE fields, one can see in their entirety 
all three zero-energy curves including the optimal path to extinction in Fig. [21(b) . Also 
shown in Fig. [2Kb) are the analytical solutions to the three zero-energy curves given by 
p — 0, q = 0, and Eq. (|20[) . There is excellent agreement between the analytical solu- 
tions of all three curves and the LCS which are found through numerical computation 
of the FTLE flow fields. 



4.2 Example 2 - SIS Epidemic Model - External Fluctuations 

As another general application of the extinction theory for finite populations, we con- 
sider the well-known problem of extinction in a Susceptible-Infectious-Susceptible (SIS) 
epidemiological model. The SIS model is given by the following system of equations: 

S =ii - nS + 7/ - (3IS, (21a) 
i = -Qi + i)I + 0IS, (21b) 

where fi represents a constant birth and death rate, /3 represents the contact rate, and 
7 denotes the rate of recovery. If we assume that the total population size is constant 
and can be normalized to S + I = 1, then Eqs. (|21a|l - (|21b[) can be rewritten as the 
following one-dimensional (ID) equation: 

i = -(» + 1 )I + pI(l-I). (22) 

The stochastic version of Eq. (|22|) is given as 

i = -(jt + i)I + 131(1 -I) + q(t) = F(I) + q (t), (23) 

where rj(t) is uncorrelated Gaussian noise with zero mean and models random migra- 
tion to and from another population ( Alonso et all [20061 : iDoering et all . 120051 ). Equa- 



tion (|23p has two equilibrium points given by I = (corresponding to the disease-free 
state) and / = 1 — (/x + 7)//? (corresponding to the endemic state). One can use the 
Euler-Lagrange equation of motion to find the optimal path of extinction from the en- 
demic state to the disease-free state, where the Lagrangian is determined by Eq. (|23ll 
and is given as follows: 

L(lJ) = [ V (t)] 2 = [i-F(I)} 2 . (24) 
Computation of the Euler-Lagrange equation gives the following: 

F(I)F'(I) -1 = 0. (25) 

If one multiplies Eq. (|25[) by I followed by an integration with respect to t, then one 
obtains 

where E is an arbitrary constant. Using the fact that the optimal path passes through 
the two equilibrium points stated above, then one finds that the optimal path to ex- 
tinction (as well as its counterpart path from the disease-free state to the endemic 
state) is given by the following equation: 



/ = ±F(I). 



(27) 
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0.2 0.4 0.6 0.8 I 
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Fig. 3 (Color online) FTLE flow fields computed using Eqs. (f28at - (l28bb with f3 = 5.0 and 
k = 1.0. The integration time is T = 5 with an integration step size of t = 0.1 and a grid 
resolution of 0.01 in both / and p. (a) Forward FTLE field with the optimal path to extinction 
given by the LCS. (b) Average of the forward and backward FTLE fields with the optimal 
path to extinction and its counterpart optimal path given by the LCS and overlaid with the 
analytical solution of the optimal paths given by Eq. 1271 . Note that the averaging affects only 
the value of the FTLE and not the structure of the FTLE field. 



As in the first example, one can numerically compute the optimal path to extinction 
using the FTLE. In this example we calculate the FTLE field using the following 2D 
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system which is equivalent to Eq. (12511 : 

7 =p, (28a) 
p =F{I)F'(I) = 09/(1 — I) — «/)(/?(! - 2/) - k), (28b) 



where k = /u + 7. Figure[3ja) shows the forward FTLE plot computed using Eqs. (|28al - 
()28bl) for T = 5, with f3 = 5.0 and k = 1.0. In Fig. one can see that the optimal 

path from the endemic state to the disease-free state is given by the ridge associated 
with the locally maximal FTLE. 

One also can find the optimal path from the disease-free state to the endemic state 
by computing the backward FTLE. By overlaying the forward and backward FTLE 
fields, one can see the optimal path to extinction along with its counterpart optimal 
path in Fig.[3jb). Also shown in Fig.[3jb) are the two analytical solutions to the optimal 
path to extinction and its counterpart optimal path which are given by Eq. Q27p . There 
is excellent agreement between the analytical solutions to the two optimal paths and 
the ridges which are found through numerical computation of the FTLE flow fields. 



4.3 Example 3 - SIS Epidemic Model - Internal Fluctuations 

We now consider the ID stochastic (internal) version of the SIS epidemic model given 
by Eq. ()22p . The probability P n (t) to observe, at time t, n infectious individuals is 
governed by the following master equation: 

Pn = {jl + 7)[(»» + 1)-Pn+1 - nP n] + £[(>» - W ~ (*» ~ - "(« - 1)3.]. (29) 



Using the formalism of lGan d 1 19871 ). one then has the following Hamiltonian associated 
with Eq. 1(25): 

H{I,p) = (/i + 7 )/(e" P - 1) + 131(1 - I)(e p - 1), (30) 
and Hamilton's equations are therefore given as: 

i=^ = -( f i + y)Ie- p + pi(l-I)e p , (31a) 
op 

V = ~ ^ = -{H + l)(e- p - 1) + /3(e p - 1)(2J - 1). (31b) 



Although there is no analytical solution for the optimal path to extinction for 
Eqs. <[5Ia ) -((3Tb ) . we can once again determine the optimal path by computing the 
FTLE flow field associated with this system. Figure [4] shows the forward FTLE plot 
computed using Eqs. (|31a[) - ((31b[) for T = 10, with /3 = 2.0 and n = 1.0. As we have seen 
previously, the optimal path to extinction from the endemic state to the disease-free 
state is given by the ridge associated with the locally maximal FTLE. 
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5 Maximal Sensitive Dependence to Initial Data Near the Optimal Path 

The main heuristic argument of this section is to show that the path which maximizes 
the probability of extinction also has a finite-time Lyapunov exponent (FTLE) that 
attains its local maximum on the path. We consider a general system of equations 
and show how the maximum unstable direction along the optimal path to extinction 
governs the local hyperbolic dynamics. 

From the Hamiltonian or Lagrangian equations of motion, the process which leads 
to extinction consists of a trajectory which emanates from a stable steady state x a £ 
ffi" and approaches the extinct state x s £ M™. Since the stable and extinct states 
are both regular saddles (or unstable foci) in the variational formulations, they both 
have hyperbolic structure. Moreover, every point along the trajectory connecting the 
two states as i ±cxj is assumed to possess a local hyperbolic structure. As an 
example, consider the Langevin problem having a vector field of position V : R" — > R n , 
which has the associated Lagrangian L(x, x) = \\x — V(x)\\ 2 /2 to describe the action. 
Converting to a Hamiltonian formulation leads one to the following 2n-dimensional 
equations of motion and Hamiltonian: 

x = p + V{x), (32a) 
p = -V'(x)p, (32b) 

H(x,p) = M-+p-V(x), (32c) 
where V'{x) = is the Jacobian matrix evaluated at x. 
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It is immediate from Eqs. (|32a[) - (|32c[) that {(x,p) \ p = 0} is an invariant manifold. 
In addition, the optimal path must lie along the H(x,p) = surface, which means that 
in addition to the p = manifold, the zero surface includes {(x,p) | p — — 2V(x)}. 

To clarify the direction along the optimal path as well as the local geometry, we 
make the following assumptions regarding V(x): 

1. V(x) is smooth, 

2. V(x a ) = V(x s ) = 0, 

3. V'(x a ) has eigenvalues with negative real parts, and V (x s ) has at least one eigen- 
value with positive real part. 

Items 2 and 3 imply that x a is an attracting steady state and x s is an unstable 
steady state in the deterministic dynamical system. We now assume that the optimal 
path must satisfy lim (x(t),p(t)) = (x s ,0), while lim (x(t),p(t)) = (a; a ,0). Since 

t— f+oo t— > — oo 

H(x(t),p(t)) — along the path, the limits provide direction along the optimal path. 
The optimal path lies on the curve 

C {x , p) ={te (-oo, oo) I p(t) = -2V(as(t))} , 

and p = corresponds to the zero fluctuation case. We shift the optimal path to the 
origin by using the following 2n-dimensional transformation: 

u = x, (33a) 
w = p + 2V(x), (33b) 

H(u, w) = J^|1L - w ■ V(u). (33c) 

The new equations of motion are now: 

u = dH/dw = w - V(u), (34a) 
w = -dH/du = V'(u)w. (34b) 



The optimal path now is described by the curve 



a 



(u,w) 



{t€(- 



3)|t»(t) = 0,u(t) = -V(u(t))}, 



while the zero fluctuation case given by p — now corresponds to w — 2V(u). 

The linearized variation along the optimal path Ci u w \ is given by the following 
matrix initial value problem from Eqs. (134all - (l34b[) : 



X = 



-V'(u(t)) In 

V'(u(t)) 



X = J(u(t),0)X, X(0) = I. 



(35) 



For a fixed time to such that u(to) — uq, the local eigenvalues of Eq. (135ft are given 
by the eigenvalues of ±V"'(«o) and are assumed to have non-zero real part for any 
(no, 0) £ C(„ w y Thus the optimal path is hyperbolic at every point. We also suppose 
the existence of a local coordinate system on the path so that there exists a set of 
linearly independent directions pointwise. 

The solution to the linear variational equation about (u, w) = (uq, 0) for 0<i<l 
is given by X(t) ~ exp (t J(uq, 0)). We assume for simplicity that the eigenvalues of 
V'(uq) have algebraic multiplicity of one. The eigenvalues of J(uq,0) are given by 
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{±\i}™ =1 , where Xi are eigenvalues of V'(uq). To examine the dynamic instability that 
dominates locally, let Amax denote the eigenvalue with largest real part and be such 
that Re(Xmax) > 0. Since —Amax has the most negative real part, its eigendirection 
denotes the strongest contracting direction. 

For any (uq, 0) on the path, the existence of a set of linearly independent eigenso- 
lutions of J(uq, 0) implies there is a transformation X = PY, with P £ L 2ri x2"^ that 
diagonalizes the linear variational equation given by Eq. (|35J1 . Without loss of gener- 
ality, we can also assume that the diagonal consists of ordered descending eigenvalues 
based on the real part. Therefore, the linear variational system has the form 



"An 



An 



For any initial value, the solution to Eq. (1361) is 



A, 



(36) 



x p (t;x ) = (xi(t),x 2 (t), ■ ■ ■ ,x 2n (t)) 



i > 
(e 



K t A 2 t 
^lo, e x 2 q,- ■■ 



, e X n Q, 
-A 2 t 



(37) 



(38) 



X {n+1)0> ' " > e 2(271-1)0' e X 2 „0). 

To show that the FTLE takes it maximum along the path, we notice that any 
point along the path is hyperbolic with a saddle structure. Therefore, we consider an 
arbitrary initial condition lying within a small domain containing the origin. Since 
almost any initial condition hits the boundary of the domain in finite time due to the 
saddle structure of the origin, we use the escape time as the final time for the FTLE. 
The definition we use of the FTLE is the direct comparison of the distance between 
two close trajectories as follows: 

cr(i;»o) = j hx{\\x p {t;xo +e) - x p (t; xq)\\), 



(39) 



where eel*. 

Defining the domain to be the 2n-dimensional hypercube D = [— 1,1]", then 
clearly any point not on the unstable manifold will escape in the x\ direction cor- 
responding to the eigenvalue with maximal real part. We exploit the fact that the 
dynamics is governed by the most unstable direction by assuming |A m ax| >> Aj|, 
i = 2 • • • n. If the initial condition lies within a distance 5 of the unstable manifold with 
< 8 « 1, then the time to escape from the domain for an arbitrary non-zero initial 
condition is given by 

;(*) 



lot! 



Amax 

Using the definition of the exponent given by Eq. ((39]), we have that 

2 i A, 1 2 I All 1 2 



(40) 



a(tf,x ) 



+ ■ 



In 



5 



+ 6~ 



: e 2 



£2n-l 



+ \St 2r . 



+ 5~ 



/ (2 hi (<5)). 



: e„\ +\S^ 



tn+1 



(41) 
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Since A ma x| >> |Aj|, and since ±A m ax dominates the expanding and contracting 
directions, then for 8 small, we may just consider 

/, „, v -A m axln(ef/<? 2 +el n 5 2 ) . , 



Furthermore, we find that 



da(t f ;x (S)) = A m ax In {e\) / \ q f 

dS 2<5(ln<5) 2 V ef 7 V ln(5 



(43) 



which can be shown to be negative assuming e\ « 1. Therefore, the FTLE as a 
function of distance to the stable invariant manifold is a decreasing function, and thus 
takes it maximum values on the manifold. 



6 Conclusions 

In this article, we have considered the dynamics of general stochastic epidemic models 
and their extinction properties in finite populations. The random fluctuations consid- 
ered were from both internal fluctuations, which arise from mass action kinetics, as 
well as external random forces, which may be due to random population migrations. 
By examining the extinction processes from the master equation perspective, eikonal 
approximations in the large population limit give a way to solve for the probability 
distribution as a function of time. A variational principle applied to the exponent of 
the probability distribution near the steady state was used to maximize the probability 
to extinction from a disease endemic state. 

Maximizing the probability to extinction means minimizing the action, which in 
turn generates a Hamiltonian (or Lagrangian) formulation that determines the flow 
from endemic to extinct states. The formulation describes the random fluctuations as 
a deterministic effective force that overcomes the instability of the extinct state. Such 
a deterministic flow describes the optimal path to extinction from an endemic state. 
Using the above variational formulation, we explicitly derived the equations of motion 
describing the optimal path to extinction in three different models from epidemiology 
as a two point boundary value problem 

The main result of our paper is that the optimal path is intimately related to the 
maximal sensitivity of the dynamics of unstable Hamiltonian flows. Specifically, we 
introduced a novel method to find the optimal path to extinction that relies on the 
computation of the finite-time Lyapunov exponent field for the dynamical flow under 
consideration. The exponents provide a measure of sensitivity to initial conditions in 
finite time. Moreover, we have shown that the system possesses maximal sensitivity near 
the optimal path to extinction. Therefore, we are able to use the finite-time Lyapunov 
exponents to dynamically evolve toward the optimal path trajectory. 

To demonstrate the equivalence of the maximal sensitivity and the optimal path 
maximizing the probability of extinction, we have considered three prototypical exam- 
ples from mathematical epidemiology. In the examples, we have considered both inter- 
nal and external noise, and we have considered both a Hamiltonian and Lagrangian 
formulation. Furthermore, in each of the three examples, we have shown that the opti- 
mal path to extinction is equated with having a (locally) maximal sensitivity to initial 
condition, which implies a relation at a fundamental level between the optimal path and 
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Fig. 5 (Color online) Snapshots taken at (a) t = 1, (b) t = 3, (c) t = 5, (d) t = 7, and (e) 
t = 9 showing the evolution of the FTLE flow field computed using Eqs. I3lat - ll3lbll with 
f) = 2.0 and n = 1.0. The integration time is T = 10 with an integration step size of t = 0.1 
and a grid resolution of 0.005 in both I and p. The final snapshot in the series is taken at 
t = 10 and is shown in Fig. [1] A video that shows the evolution of the FTLE flow field can be 
found online as ancillary material. 



the FTLE. Even though there exist many possible paths to extinction, the dynamical 
systems approach converges to the path that maximizes the probability to extinction. 

An example of the evolution of the FTLE flow field showing the convergence of the 
locally maximal FTLE ridge to the optimal path is shown in Fig. [5] The FTLE flow 
field is computed for T = 10 using Eqns. (|31a[) - (l31b|l for the SIS epidemic model with 
internal fluctuations (Sec. 14.3)) . Figure [S] shows snapshots of the FTLE field taken at 
t = 1, t = 3, t = 5, t = 7, and t = 9. The final snapshot at t = 10 is shown in Fig. [4] A 
video that shows the evolution of the FTLE flow field can be found online as ancillary 
material. 

The parameter values chosen for the three examples are such that the extinct and 
endemic states are far away from one another, implying that the system is operating 
far from any bifurcation points. This result is very importa nt, since in gen e ral, n o 
approximate analytical treatment, like the one performed in iDvkman et alj ( 20081 ) . 



is possible if the system's dynamics is not sufficiently close to the bifurcation point. 
However, any scaling behavior of the exponent in the probability of extinction may 
still be computed along the optimal path, which is the important advance afforded by 
the new procedure proposed in this paper. 

In the future, we plan on considering more com plicated systems. Of parti cular 
interest are bistable systems (e .g. adaptive network s (|Shaw and Schwartl I2008T ) and 
the Schlogi birth-death process dDoering et all 120071) '). and higher-dimensional systems 



(e.g. multistrain epidemic models~ (|siiawet'ld] . I2007I ) ) . Because the method is general, 
and unifies dynamical systems theory with the probability of extinction, we expect that 
any system found in other fields can be understood using this approach. 
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